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Summary. Combining deeper insight of Einstein's equations with sophisticated 
numerical techniques promises the ability to construct accurate numerical imple- 
mentations of these equations. We illustrate this in two examples, the numerical 
evolution of "bubble" and single black hole spacetimes. The former is chosen to 
demonstrate how accurate numerical solutions can answer open questions and even 
reveal unexpected phenomena. The latter illustrates some of the difficulties en- 
countered in three-dimensional black hole simulations, and presents some possible 
remedies. 



1 Introduction 

Extracting the full physical content from Einstein's equations has proven to 
be a difficult task. The complexity of these equations has allowed researchers 
only a peek into the rich phenomenology of the theory by assuming special 
symmetries and reductions. Computational methods, however, are opening 
a new window into the theory. To realize the full utility of computational 
solutions in exploring Einstein's equations, several questions must first be 
addressed. Namely, a deeper understanding of the system of equations and 
its boundary conditions, the development and use of more refined numerical 
techniques and an efficient use of the available computational resources. 

In recent years, considerable advances have been made in some of these 
issues, allowing for the analysis of complex physical systems which arguably 
must be tackled numerically. In the present article we highlight some recent 
analytical and numerical techniques and apply them to two practical appli- 
cations. The first application is the numerical evolution of bubble spacetimes 
in five-dimensional Kaluza-Klein theory. We study their dynamical behavior, 
the validity of cosmic censorship -in a set-up which a-priori would appear 
promising to give rise to violations of the conjecture- and reveal the exis- 
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tence of critical phenomena. As a second application, we discuss the numeri- 
cal evolution of single black hole spacetimes. Here we consider some analytical 
and numerical difficulties in modeling these systems accurately. We discuss a 
method to alleviate some of these problems, and present tests to demonstrate 
the promise of this method. 

2 Analytical and Numerical tools 

In the Cauchy formulation of General Relativity, Einstein's field equations are 
split into evolution and constraint equations. Numerical solutions are found 
by specifying data on an initial spacelike slice, subject to the constraints, and 
by integrating the evolution equations to obtain the future development of 
the data. Owing to finite computer resources, one is forced to use finite, and, 
in practice, rather small computational domains to discretize the problem. 
This raises several important issues. 

The fundamental property for any useful numerical solution is that the 
solution must convergence to the continuum solution in the limit of infinite 
resolution. A prerequisite for a well-behaved numerical solution is a well- 
posed continuum formulation of the initial-boundary value problem. In cer- 
tain cases, the well-posed continuum problem can then be used to construct 
stable numerical discretizations for which one can a priori guarantee conver- 
gence. In particular, this can be achieved for linear, first-order, symmetric 
hyperbolic systems with maximally dissipative boundary conditions ["QI21E]- 
This is briefly discussed in Sec. 12. II for a detailed description and an extension 
to numerical relativity see Refs. 0JE1E1C1- 

The application of these ideas in general relativity is, naturally, more com- 
plicated. First, Einstein's equations are nonlinear and so it is much harder 
to a priori prove convergence. However, a discretization that guarantees sta- 
bility for the linearized equations should already be useful for the nonlinear 
equations, especially for those systems with smooth solutions as expected for 
the Einstein equations when written appropriately. This is because in a small 
enough neighborhood of any given spacelike slice, the numerical solution can 
be modeled as a small amplitude perturbation of the continuum solution. 

The constraint equations in general relativity bring additional compli- 
cations and greatly restrict the freedom in specifying boundary and initial 
data. This is illustrated and further discussed in Section l2~2l Section [2*31 dis- 
cusses issues regarding the stability of the constraint manifold. The manifold 
is invariant with respect to the flow defined by the evolution system in the 
continuum problem. Numerically, however, small errors in the solution aris- 
ing from truncation or roundoff error may lead to large constraint violations 
if the constraint manifold is unstable. Section 12.31 discusses a method for 
suppressing such rapid constraint violations. 
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2.1 Guidelines for a stable numerical implementation. 

A simple numerical algorithm, or "recipe," can be followed to solve first 
order, linear symmetric hyperbolic equations with variable coefficients and 
maximally dissipative boundary conditions, for which stability can be guar- 
anteed. It is based on finite difference approximations with spatial difference 
operators that satisfy the summation by parts (SBP) property. This property 
is a discrete analogous of integration by parts, which is used in the derivation 
of energy estimates, a key ingredient for obtaining a well posed formulation 
of the continuum problem. SBP allows to obtain similar energy estimates for 
the discrete problem. 

Employ spatial difference operators that satisfy SBP on the computational 
domain. For the sake of simplicity, consider a set of linear, first order sym- 
metric hyperbolic equations in the one-dimensional domain x € (a, b) which 
is discretized with points Xj — a + jAx, j — . . . N, where Ax = (b — a) /N. 
Now let us introduce the discrete scalar product, 

N 

(u,v) := Ax 2J VijUiVj , (1) 

i,3=0 

for some positive definite matrix with elements <ry which in the contin- 
uum limit Ax — y approaches the Li norm (u, v) := j^uvdx. At the 
continuum level, the derivative operator d/dx and scalar product satisfy 
integration by parts, i.e. (du/dx,v) + (u,dv/dx) = uvy a , which in the dis- 
crete case is translated into a finite difference operator D which satisfies 
(Du,v) + (u,Dv) = uv\ b a and approaches d/dx in the continuum limit. The 
simplest difference operator and scalar product satisfying SBP are 

Du = (uj+i — Uij/Ax , coo — h for i = 

Du = (u i+1 - Ui-{)/{2Ax) , era = 1 for i = 1 . . . N - 1 (2) 

Du = (ui — Ui-i)/Ax , <?nn = \ for i — N 

where the scalar product is diagonal: cry = for i ^ j. Higher order op- 
erators satisfying SBP have been constructed by Strand Additionally, 
when dealing with non-trivial domains containing inner boundaries, addi- 
tional complexities must be addressed to attain SBP, see Ref. @]. The finite 
operator D is then used for the discretization of the spatial derivatives in the 
evolution equations, thus obtaining a semi-discrete system. 

Impose boundary conditions via orthogonal projections "JSL This ensures the 
consistent treatment of the boundaries, guaranteeing the correct handling 
of modes propagating towards, from and tangential to the boundaries. An 
energy estimate can be obtained for the semi-discrete system. 

Implement an appropriate time integration algorithm. The resulting semi- 
discrete system constitutes a large system of ODE's which can be numerically 
solved by using a time integrator that satisfies an energy estimate [HI E] ■ 
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Consider adding explicit dissipation It is well known that finite difference 
approximations do not adequately represent the highest frequency modes 
on a given grid, corresponding to the shortest possible wavelengths that 
can be represented on the grid. If the smallest spacing between points is 
A, the shortest wavelength is A m , n = A with the corresponding frequency 
kmax = ^Tr/^min- These modes can, and often do, travel in the wrong direc- 
tion. For this reason, it is sometimes useful to add explicit numerical dissipa- 
tion to rid the simulation of these modes in a way that is consistent with the 
continuum equation at hand. As finer grids are used, the effect of this dissi- 
pation becomes smaller and acts only on increasingly higher frequencies. The 
dissipation operators are constructed such that discrete energy estimates, ob- 
tained using SBP, are not spoiled. Explicit expressions for such dissipation 
operators are presented in Ref. @]. 

To summarize, beginning with a well-posed initial-boundary value prob- 
lem, we mimic the derivation of continuum energy estimates for the discrete 
problem using (1) spatial derivative operators satisfying summation by parts, 
(2) orthogonal projections to represent boundary conditions and (3) choosing 
an appropriate time integrator. 

2.2 Constraint-preserving boundary conditions 

As discussed above, a numerical implementation of any system of partial dif- 
ferential equations necessarily involves boundaries. Unless periodic boundary 
conditions can be imposed, as is often the case for the evolution on compact 
domains without boundaries, one deals with an initial-boundary value prob- 
lem, and thus has to face the question of how to specify boundary conditions. 
In theories that give rise to constraints, like general relativity, such conditions 
must be chosen carefully to ensure that the constraints propagate. 

As a very simple illustration, consider the Id wave equation u jtt = u^ xx 
on the half line x > 0. Let us reduce it to first order form by introducing the 
variables / = u x and g = u t — bu^ x , with b a negative constant: 



At the boundary x = 0, the system has two ingoing fields, given by u and 
Vi n = g + bf — f, and one outgoing field. However, the ingoing fields cannot 
be given independently, as we see next. The constraint C = f — u tX = 
propagates as C.t = bC tX and so C is an ingoing field with respect to x = 0. 
Therefore, we have to impose the boundary condition C = which implies 
the condition u x = / on the main variables. We can replace this with a 
condition that is intrinsic to the boundary by using the evolution equation 
||2J in order to eliminate the ir-derivative and obtain 



u tt = bu, x + g, 

g,t = -6 5 ,x + (i-6 2 )/,. 

f,t = 9,x + bf, x . 



(3) 
(4) 
(5) 
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u t = bf+g, (6) 

This equation provides an evolution equation for determining u at the bound- 
ary, which guarantees that the constraint C = is preserved throughout 
evolution. It can be complemented by the Sommerfeld condition Vi n = 0. 

This simple example gives just a glimpse of the different issues involved 
in prescribing constraint-preserving boundary conditions. The case of Ein- 
steins's field equations is more complicated; we refer the interested reader to 
Refs. [IlEIlinillSIHnSlEinilllllllllllinilll A major difficulty is 
the fact that, in general, constraint-preserving boundary conditions do not 
have the form of maximal dissipative boundary conditions, and for this rea- 
son it has proven to be difficult to find well posed initial-boundary value 
formulations of Einstein's equations that preserve the constraints. 

2.3 Dealing with "too many" formulations. Parameters via 
constraint monitoring 

Formulations of the Einstein equations are often cast in symmetric hyperbolic 
form by adding constraints to the evolution equations multiplied by parame- 
ters or spacetime functions. The symmetric hyperbolicity condition partially 
restricts these parameters, however, considerable freedom in the formulation 
exists in choosing these free parameters (see, for instance, |28|). Analytically, 
when data are on the constraint surface, all allowed values for these parame- 
ters are equally valid. Off of the constraint surface, however, different values 
of these parameters may be regarded as representing "different" theories. It is 
no surprise then that numerical simulations are sensitive to the values chosen 
for these parameters, as numerical data rarely are on the constraint surface. 
Unfortunately, the parameters in current simulations are proving to be ex- 
tremely sensitive. Relatively mall variations in these parameters (within the 
allowed range for a symmetric hyperbolic formulation) produce run times in 
simulations that vary over several orders of magnitude, as measured by an 
asymptotic observer. 

Furthermore, the parameters are not unique. Values convenient for one 
physical problem might be inappropriate in another. Recently, a method to 
dynamically choose these parameters — promoted to functions of time — was 
introduced that naturally adapts to the physical problem under study |2*3|. 
Basically, one exploits the freedom in choosing these functions to control the 
growth rate of an energy norm for constraint violations. Since this norm is 
exactly zero analytically, this provides a guide to choosing the parameters 
that will drive the solution to one that satisfies the constraints. This method 
provides a practical solution to this problem of choosing parameters, although 
it may not be the most elegant solution. Ideally, one would like to understand 
how the growth rate of the solution depends on the these parameter values 
in order to choose them appropriately. This would require sharp growth es- 
timates, however, which are still unavailable. While further understanding is 
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gained in this front, this practical remedy can be of much help in present 
simulations. We summarize here the essential ideas of this method. 

Consider a system of hyperbolic equations with constraint terms, C c , writ- 
ten schematically as 



i a = '^2 ^ 6 ( u ' *) x )dbU a + B a (u, t, x) + fi ac C c (u, dju) , 



(7) 



where u a , B a and C' c are vector valued functions, and /i ac is a matrix (gen- 
erally not square) that is a function of the spacetime (C c represents a vector 
function of general constraint variables). The indices {a, 6, c} range over each 
clement of the vector or matrix functions, while the indices {i,j,k} label 
points on a discrete grid. We define an energy or norm of the discrete con- 
straint variables as 



Af(t) 



1 



e 2 l fl X TlyTl Z 



EE^w 2 



(8) 



where n x , n y , n z are the number of points in each direction. The grid indices 
k} are suppressed to simplify the notation. The time derivative of the 
norm can be calculated using Eq. Q 

= jhom + Tr ( M jM) ) (9) 
and therefore can be known in closed form provided the matrix valued sums 



•J-tlOTTl 



c„ 



T'' 



EE 

ij k a,b 

c 

c a 



du b 



B h 



^ dD k u a 



(10) 



EE 

ijk a 
dCg 

du b 



V- 9C a 



t-r' dD k u b 

k 



Cc 



(11) 



are computed during evolution. Here Di is the discrete derivative approxi- 
mation to di. We then use the dependence of the energy growth on the free 
constraint-functions to achieve some desired behavior for the constraints, i.e., 
solving Eq. © for fi ac . For example, if we choose 6 



Jif = -aAf, 



a > 0, 



(12) 



There is a slight abuse of notation here, in the sense that a does not denote 
an index, as before. Similarly, the subscript in n a indicates that the quantity is 
related to a through Eq. 11411 . 
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any violation of the constraints will decay exponentially 

Af(t + At)=N(t)e- aAt . (13) 

As discussed in Ref. one good option among many others seems to be 
choosing a tolerance value, T, for the norm of the constraints that is close to 
the initial discrete value, and solving for [i ac such that the constraints decay 
to this tolerance value after a given relaxation time. This can be done by 
adopting an a such that after some time r = n a At the constraints have the 
value T. Replacing Af(t + At) by T in equation and solving for a gives 

If one then solves 

Af = -aAf = I hom + traced x 1^) (15) 

for n, with a given by Eq. I|14|l . the value of the norm Af(t + r) should be 
T, independent of its initial value. Therefore, Eq. (|T5|) serves as a guide to 
formulate a practical method to choose free parameters in the equations with 
which the numerical solution behaves well with respect to the satisfaction 
of the constraints. Naturally, if one deals, as it is often the case, with more 
than one free parameter, Eq. (|15fl must be augmented with other conditions 
to yield a unique solution. This extra freedom is actually very useful in pre- 
venting large time-variations in the parameters that are sometimes needed 
in order to keep the constraints under control. These large variations do not 
represent a fundamental problem but a practical one, due to the small time 
stepping that they require in order to keep errors due to time integration 
reasonably small. One way to prevent this is by using this extra freedom to 
pick up the point in parameter space that not only gives the desired con- 
straint growth, but also minimizes the change of parameters between two 
consecutive timesteps. 

Rather than including the full details on the particular way we have im- 
plemented the method, we describe here a simple example to illustrate its 
application. Assume, for instance, that within a particular formulation only 
two free functions, {k,uj}, are employed, Eq. ijTB|) formally evaluates to 

AY= -aAf = X hom + kZ k + ul" . (16) 

Now, we exploit the freedom in the free functions to adjust the rate of change 
of the energy Af if the values of \I hom ) 2 K , jy } are known. In practice, these 
are easily obtained during evolution. Once these are known, Eq. (|16f) . coupled 
to the requirement that vary as little as possible from one evaluation 

to another, results in a straightforward strategy to evaluate preferred values 
of the free parameters. This is done at a single resolution "test" run and, 
through interpolation in time, continuum, a priori defined parameters which 
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keep the constraints under control for the given problem are obtained. De- 
pending on the formulation of the equations, the free parameters might have 
to satisfy some conditions in order for symmetric hyperbolicity to hold, which 
can restrict the range of values these parameters can take. Nevertheless, even 
within a restricted window, the technique allows one to adopt the most con- 
venient values these parameters should have for the problem at hand. 

3 Applications 

We now present applications of the techniques previously discussed. The goal 
is to illustrate how well-resolved simulations can indeed serve as a powerful 
tool to understand particular problems. To this end we have chosen a problem 
found in higher dimensional general relativity. A second application is that of 
the simulation of single black hole spacetimes, where the issue of the a priori 
lack of a preferred formulation is illustrated. 

3.1 Bubble spacetimes 

As a hrst application we concentrate on the study of bubble spacetimes and 
elucidate the dynamical behavior of configurations with both positive and 
negative masses and their possible connection to naked singularities. Bubble 
spacetimes have been studied extensively within five-dimensional Kaluza- 
Klein theory. These are five-dimensional spacetimes in which the circum- 
ference of the "extra" dimensions shrinks to zero on some compact surface 
referred to as the "bubble". These bubbles were initially studied by their 
relevance in the quantum instability of flat spacetime |25j . as bubbles can 
be obtained via semi-classical tunneling from it. They were later extended to 
include data corresponding to negative energy configurations (at a moment 
of time symmetry) [261127) . As mentioned, among the reasons for considering 
negative energy solutions is that naked singularities are associated with them. 
Therefore, these solutions are attractive tests of the cosmic censorship conjec- 
ture. Additionally, bubble spacetimes can also be obtained by double- Wick 
rotation of black strings, whose stability properties (or lack thereof) have been 
the subject of intense scrutiny in recent years. These features make bubble 
spacetimes both interesting an relevant for gravity beyond four-dimensions, 
and thus attention has been devoted to fully understand their behavior. As we 
will see, even when the "analytical" study of the problem is greatly simplified 
by symmetry assumptions, many lingering questions remain and numerical 
simulations provided a viable way to shed light into them. Furthermore, these 
simulations were also key to 'digging out' a few unexpected features of the 
solution. 

In order to obtain a complete description of the dynamical behavior of 
these spacetimes, a numerical code, implementing Einstein equations in 5D 
settings, and capable of handling the possibly strong curvature associated 
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need be constructed. Fortunately, the assumption of a SO(3) x U(l) sym- 
metry simplifies the treatment of the problem, which can be reduced to a 
1 + 1 manifold. This, in turn, renders the problem quite tractable by the cur- 
rently available computational resources, though as we will see, considerable 
care must be placed at both analytical and numerical levels for an accurate 
treatment of the problem. 

Initial data 

We consider a generalization of the time symmetric family of initial data 
presented in We start with a spacetime endowed with the metric 

dr 2 

ds 2 = -dt 2 + U(r)dz 2 + — — — + rW, (17) 

U(r) 

where dfl 2 = dd 2 +sm 2 fldip 2 is the standard metric on the unit two-sphere S 2 
and U (r) is a smooth function that has a regular root at some r = r + > 0, 
is everywhere positive for r > r + and converges to one as r — > oo. The 
coordinate z parameterizes the extra dimension S 1 which has the period 
P = 47r/f7'(r+). The resulting spacetime {t,z,r > r+,#, <p\ constitutes a 
regular manifold with the topology Rx R 2 x S 2 . The bubble is located where 
the circumference of the extra dimension shrinks to zero, that is, at r = r + . 

Additionally, we consider the presence of an electromagnetic field of the 
form 

- dx» A dx v = dy(r) A dz, (18) 

where j(r) is a smooth function of r that converges to zero as r — > oo. 
The symmetries of the problem would also allow for a non-trivial electric 
component of the field. However, it is not difficult to show that Maxwell's 
equations imply that such a field necessarily diverges at the location of the 
bubble. For this reason, in the following, we only consider the case of vanishing 
electric field. 

In this article, we consider initial data with 

>T(r) = k{r+ n -r- n ), (19) 

where k is an arbitrary constant and n an integer greater than one. This 
field generalizes the ansatz considered in |57], where only the case n = 2 
was discussed, and allows for different interesting initial configurations. In 
the time-symmetric case, initial data satisfying the Hamiltonian constraint 
obeys 

m b k 2 , 

u ( r ) = 1 - 7 + ^ - . ( 2 °) 

with k = kn/y/(n — l)(2n — 1) and free integration constants m and b. Here, 
the parameter m is related to the ADM mass via Madm = m/4. The fact 
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that the bubble be located at r = r + requires that = U(r + ) = l — m+b — k 2 

I, k = k/rl . 

< r+U'(r+) = 2 - m + 2(n - l)k 2 (21) 



where m = m/r + , b = b/r 2 _, k = fc/r" . We also require 



and avoid the conical singularity at r = r + by fixing the period of z to 
P = An /U 1 (r+). It can be shown that the initial acceleration of the bubble 
area A with respect to proper time is given by 



A = 8ir 



Ak 2 

1 - to — (n - l)(n - 2) 

o 



(22) 



For n = 2, as discussed in Ref. |28| . this implies that negative mass bubbles 
start out expanding (the initial velocity of the area is zero since we only con- 
sider time-symmetric initial data), while for large enough positive mass the 
bubble starts out collapsing. In the vacuum case, our numerical simulations 
suggest that initially collapsing bubbles undergo complete collapse and form 
a black string. In the non-vacuum case however, the strength of the elec- 
tromagnetic field can modify this behavior completely. We will see that for 
small enough k the bubble continues to collapse whereas when k is large the 
bubble area bounces back and expands. Interesting behavior is obtained at 
the critical value for k which divides the phase space between collapsing and 
expanding solutions. 

For n > 2 it is possible to obtain initial configurations with negative mass 
and negative initial acceleration (2H|. This can potentially give rise to a col- 
lapsing bubble of negative energy, and thus to a naked singularity. However, 
our numerical results (2Hj suggest that cosmic censorship is valid: The bubble 
bounces back and starts out expanding. 



Equations 

In order to study the time evolution of the initial data sets given on a 
t = const slices of the metric ()17|l and the electromagnetic field (I18JI . it 
is convenient to introduce a new radial coordinate R — R(r) which facilitates 
the specification of regularity conditions at the bubble location. This new 
coordinate is defined by 

R(r) = yV 2 - r\ , r>r+ . (23) 

The metric (|17|l now reads 

ds 2 = ~a 2 dt 2 + e 2a dR 2 + 2 R \ 9 e 2b dz 2 + (r 2 + + R 2 )e 2c dH 2 , (24) 

7 + + R z 

with a = 1, e' 2a = e 2b = {r\ + R 2 )U{R)/R 2 , c = 0. Since U(R) = const ■ 
(R/r+) 2 + 0{R i ) near R — 0, and U(R) converges to one in the asymptotic 
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region, a and b are regular functions. An explicit example is the initial data 
corresponding to the zero mass Witten bubble where U = 1 — (r+/r) 2 
and thus a = b = 0. When studying the time evolution of the initial data 
sets discussed above, we consider the metric <|24[) where a, a, b and c are 
functions of t and R. As we will see, the coordinate R is well suited for 
imposing regularity conditions at the bubble location since (i?, z) represent 
polar coordinates near the bubble, R = being the center, and z assuming 
the role of the angular coordinate. In order to avoid a conical singularity, z 
must have the period 2-7rr + e a ~ b . For this to be constant we need to impose 
the boundary condition a(t, 0) — b(t, 0) = const at R = 0. 

Similarly, the electromagnetic field l|18|) is written in the form 

1 R 

- F^v dx 11 A dx v = e b {■K 1 dt + d 7 dR) A dz, (25) 



2 



R 2 



where the functions 7r 7 and d 7 depend on t and R and satisfy 7r 7 = and 
eL = e~ b <9 r 7 at the initial time. 

We choose the following gauge condition for the lapse 

log(a) = a + X(b + 2c), (26) 

with a parameter A which, in our simulations, is either zero or one. For A = 1 
the resulting gauge condition is strongly related to the densitized lapse con- 
dition often encountered in hyperbolic formulations of Einstein's equations: 
Indeed, the square root of the determinant of the four metric belonging to 



Eq. <|22Jl is given by VgW = e a+b+2c R^J r\ + R 2 sintf, so Eq. (gSJl sets a 
equal to the square root of the determinant of the four metric but divides the 



result by the factor Ry r \ + R 2 sin$ which is singular at the bubble, at the 
poles d = 0, 7r/2 and in the asymptotic region. For A = 0, the condition H26fl 
implies that the two-metric —a 2 dt 2 + e 2a dR 2 is in the conformal flat gauge. 
As we will see, the principal part of the evolution equations is governed by 
the d'Alembertian with respect to this metric. Since the two-dimensional 
d'Alembertian operator is conformally covariant, the resulting equations are 
semi-linear in that case. In particular, this implies that the characteristic 
speeds do not depend on the solution that is being evolved. 

The field equations resulting from the five-dimensional Einstein-Maxwell 
equations split into a set of evolution equations and a set of constraints. The 
evolution equations can be written as 

A = e- XF \{A' + 2G')e x( - 2B+F A' - 3(A - 1)(C + G') 2 e 2XB - (A + l)V 
+ 2XAB - A(A + 1)B 2 - 3(A + 1)C 2 + G [(1 - X)n 2 - (1 + A)e 2AS 4^7) 



g _ e (A-l)B-(A+2)F 



B i e (\+1)B+(X+2)F 



ir 2 + + 2R 2 
{r 2 + +R 2 ) 2 
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+ (A - l)B 2 - 2V, 
C = e (A-l)B-F h C > + G ') e (A+l)B+F' 

„2AB j21 



V+(X-1)BC, 



2.G r 9 

_ AB-2(C+G) 



r2 +i?2 

1? — ' 



d ^XB+2(G+G) 



-(B-2C) 



R 



(XB - 2C)vr 7 , 

„B-2C 



R 2 



(B - 2C)d 1 



(28) 

(29) 
(30) 

(31) 



where we have set A 



Xb + 2(A + l)c, B = b + 2c, C = c and G 



log(rl + R 2 )/2, F = log(i?) + G and V = e 2 ^" 30 '/^ + R 2 ). Here, and in 
the following, a dot and a prime denote differentiation with respect to t and 
i?, respectively. The evolution equations constitute a hyperbolic system on 
the domain R > 0. 

The constraints are the Hamiltonian and the R component of the mo- 
mentum constraint, given by C — 0, Cr — 0, where 



c = e (A-l)B-(A+2)F 

' 3r 2 + 2R 2 



D (\+l)B+(\+2)F 



R 2 ) 2 



(B' + F'){A' + 2G') + 3(C" + G') 2 



„2AB 



V - (A - XB)B + 3C 2 + G [it 2 + e 2XB d 2 



r — P A-2C 

— e 

+ 2Gir~ f d~ f . 



(32) 

- (B' + F') \A - (A + 1)B] + 2(C + G'){3C - B) 

(33) 



Regularity conditions 



The evolution equations contain terms proportional to e~ F which diverge like 
1/R near R = 0, and therefore, regularity conditions have to be imposed at 
R = 0. This is achieved by demanding the boundary conditions 



A' = B' = C' = it-, = at R = 0. 



(34) 



Assuming that the fields are smooth enough near R — 0, it then follows that 
the right-hand side of the evolution equations is bounded for R — > 0. Next, 
as discussed above, the avoidance of a conical singularity at R — requires 
that A — (A + 1)B = a — b must be constant at R = 0. We show that this 
condition is a consequence of the evolution and constraint equations, and 
of the regularity conditions Ij34(l . Using the evolution equations in the limit 
R — > and taking into account the conditions (|34J) . we find 



0, 



— A)B 



(A + l)sH = -(X + l)e^-^ B C 

i J R=0 



R=Q 



(35) 
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This means that if the Hamiltonian constraint is satisfied at R — (or in 
the case that A = — 1 even if the constraints are violated), the condition 
A — (A + 1)B R=0 = const will hold provided that the initial data satisfies 

A — (A + 1)B = 0. Next, we analyze the propagation of the constraint 
R=o 

variables C and Cr and show that the regularity conditions (|34[1 and the evo- 
lution equations imply that the constraints are satisfied at each time provided 
they are satisfied initially. 



Propagation of the constraints 

First, we notice that the vanishing of the momentum constraint requires that 
A — (\ + 1)B =0 because of the factor F' which diverges like 1/R near 

R=0 

R = in the definition of Cr. This is precisely the condition a(t, 0) — b(t, 0) = 
const discussed above. However, for this condition to hold, we first have to 
show that the momentum constraint actually vanishes. In order to see this, 
we regularize the constraint variables and define C — e F C, Cr — c f Cr. Now 
the regularity conditions (|34|l imply that Cr is regular and that C vanishes at 
R = 0. As a consequence of the evolution equations and Bianchi's identities, 
the constraint variables obey the following evolution system 



d C = e Wd R e^ B C 



+ (3A - 1)BC, (36) 



8 t C R = e-^ B - XF d R \e^ B+XF c] + (A - l)BC 



(37) 



which is regular at R — 0. Defining the energy norm 



m = \ 



; 2S+AFJ2 + e 2{\ + l)B+\F£2\ dR 



(38) 



taking a time derivative and using the equations l|36|l . (|37|l we obtain 



dt 



._>,a-i,/.m-\/ Tc , ; , ' -.. \ I B(3e 2B+XF C 2 + 2e 2 ^ B + XF C 2 R 'jdR. 

(39) 

The boundary term vanishes because of the regularity conditions at R = 
and under the assumptions that all fields fall off sufficiently fast as R — » oo. 
If B is smooth and bounded, we can estimate the integral on the right-hand 
side by a constant C times £, and it follows that £{t) < e ct £(0). This shows 
that if the constraints are satisfied initially, they are also satisfied for all t > 
for which a smooth solution exists. In the gauge where A = we even obtain 
the result that the norm of the constraints cannot grow in time. 

To summarize, the boundary conditions (|34|l imply that the constraints 



C = 0, 
tion. 



Cr = and A - (A + l)B 



R=0 



= are preserved throughout evolu- 
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Outer boundary conditions 

For numerical computations, our domain extends from R = to R = R ma x 
for some Rmax > 0. Now we have to replace the estimate l|39|) by the estimate 



d 
dt 



£ = e 2(A+l)B+AF+Ao C ~ Cfl 



+ C£, 



(40) 



and it only follows that the constraints are zero if we control the boundary 
term at R = R m ax- For this reason, we impose the momentum constraint, 
Cr = 0, at R = R m ax- This condition results in an evolution equation for B' 
at the outer boundary. We combine this condition with the Sommcrfcld-likc 
conditions at R = Rmax, 



A + A' = 0, C + C = 0, tt 7 + (L 



0. 



(41) 



Numerical implementation 

Next, we discuss the numerical implementation of the above constrained evo- 
lution system. In order to apply the discretization techniques discussed in 
Sect. El we first recast the evolution equations into first order symmetric hy- 
perbolic form by introducing the new variables tta = A, ttb = B, ttc = C 
and dA = A' + 2G" , d B = B' , d c = C + G' . The resulting first order system 
is then discretized by the method of lines. Let us first discuss the spatial 
discretization which requires special care at R = because of the coefficients 
proportional to 1/R that appear in the evolution equations. To this end, 
consider the following family of toy models 

TV = R 1 ~ n dii(R n ~ 1 d), (42) 
d = d R TT, (43) 

where R > is the radial coordinate, and n = 1,2,3,... We impose the 
regularity condition d — at R — 0, which, for sufficiently smooth fields, 
implies that 7r = ndud at R = 0, and assume that the fields vanish for R 
sufficiently large. The toy model (|42II43|I corresponds to the n-dimensional 
wave equation for spherically symmetric solutions. The principal part of our 
evolution system has precisely this form near R = 0, where n is given by A + 1, 
A + 3, 2, 1 for the evolution equations for tta, t^b, t^c and 7r 7 , respectively. 
The system l|42H43|l admits the conserved energy 

1 f°° 

E=~ R n - 1 (tt 2 + d 2 ) dR. (44) 

A second order accurate and stable numerical discretization of the system 
can be obtained as follows: We assume a uniform grid Rj — jAR, 
j = 0, 1, 2..., approximate the fields tx and d by grid functions 77j — tt(R = 
Rj), dj — d(R — Rj), and consider the semi-discrete system 
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ttj = R)~ n D {R n - 1 d) for j > and tt = — di , (45) 

J laR 

dj = D TTj for j > and d = , (46) 

where for a grid function Uj, (Z?o u )j = i u j+i ~ Uj-i)/(2AR) is the second 
order accurate centered differencing operator. It is not difficult to check that 
this scheme preserves the discrete energy 

E discrete = ^ £ R V (*i + 4) + ^ R r^o (47) 
i=i n 

which proves the numerical stability of the semi-discrete system. Finally, we 
use a third order Runge-Kutta algorithm in order to perform the integration 
in time. By a theorem of Levermore [5], this guarantees the numerical stability 
of the fully discrete system for small enough Courant factor. 

We apply these techniques for the discretization of our coupled system. 
The outer boundary conditions are implemented by a projection method. Of 
course, the resulting system is much more complicated than the simple toy 
model problem presented above, and we have no a priori proof of numerical 
stability. Nevertheless, we find the above analysis useful as a guide for con- 
structing the discretization. Our resulting code is tested by running several 
convergence tests, and its accuracy is tested by monitoring the constraint 
variables C and Cr and the quantity A — (A + 1)B = 0. 



Results 

Here we discuss the results for the numerical evolution of the initial data 
defined by Eqs. (|17I - I20|) . We start by reviewing the evolution of the initially 
expanding bubbles and the initially collapsing negative mass bubbles and 
then focus on the initially collapsing positive mass bubbles. 

Brill- Horowitz initially expanding case 

The Brill-Horowitz initial data (n = 2) in the case of vanishing electromag- 
netic field is evolved. The bubble area A as a function of the proper time r 
at the bubble is shown in Fig. ^ for different values of the mass parameter 
m. As expected, the lower the mass of the initial configuration, the faster the 
expansion. Empirically, and for the parameter ranges used in our runs, we 
found that at late times the expansion rate obeys 

A 2-fh ,„ . 

A"7&=oy (48) 

where a dot denotes the derivative with respect to proper time r. In particular 
this approximation is valid for the bubble solution exhibited by Witten [2"5] 
which describes the time evolution in the case m = 0. 
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Fig. 1. Bubble area vs. proper time at the bubble. In this and the following plots, 
we set r+ = 1. The figure shows four illustrative examples of bubbles whose initial 
acceleration is positive. As it is evident, the expansion of the bubble continues and 
the difference is the rate of the exponential expansion. The relative error in these 
curves, estimated from the appropriate Richardson extrapolated solution in the 
limit A -> 0, is well below 0.001%. 

Collapsing negative mass case 

We here restrict to cases with negative masses that start out collapsing. In- 
terestingly enough we find that even when starting with large initial negative 
accelerations, which in turn make the bubble shrink in size to very tiny val- 
ues, it bounces back without ever collapsing into a naked singularity. As an 
example, Fig. [21 shows the bubble's area versus time for different values of n 
and k. The initially collapsing bubbles decrease in size in a noticeable way but 
this trend is halted and the bubbles bounce back and expand. Although we 
have not found a simple law as that in Eq. (|48|> . clearly the bubbles expand 
exponentially fast. Therefore, it seems not to be possible to "destroy" the 
bubble and create a naked singularity. This situation is somewhat similar to 
the scenarios where one tries to "destroy" an extremal Reissner-Nordstrom 
black hole by attempting to drop into it a test particle with high charge 
to mass ratio. There, the electrostatic repulsion prevents the particle from 
entering the hole |3l)|. 

Brill- Horowitz initially collapsing case 

Next, we analyze the Brill-Horowitz initial data for the case in which the 
bubble is initially collapsing (notice that for n = 2 this implies that the 
ADM mass is positive). While our numerical simulations reveal that in the 
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Fig. 2. Bubble area vs. proper time at the bubble. The figure shows three il- 
lustrative examples of bubble with negative mass (m = —0.1 each) whose initial 
acceleration is negative. As it is evident, the collapse of the bubble is halted and 
the trend is completely reversed. The error in these curves is estimated to be well 
below 0.001%. 

absence of the gauge field such a bubble continues to collapse, we also show 
that when the gauge field is strong enough, the bubble shrinks at a rate which 
decreases with time and then bounces back. 

Obviously, if the collapse trend were not halted, a singularity should form 
at the origin. Since the ADM mass is positive, one expects this singularity to 
be hidden behind an event horizon, and one should obtain a black string. In 
fact, for the solutions which are initially collapsing and which have vanishing 
gauge field, we observe the formation of an apparent horizon. Furthermore, 
we compute the curvature invariant quantity Ir AH at the apparent horizon 
(as discussed in where I = R a b c dR abcd is the Kretschmann invariant and 
tah the areal radius of the horizon. For a neutral black string, this invariant 
is 12. Figure El shows how this value is attained after the apparent horizon 
forms for representative vacuum cases (with m = 1.1 and m = 1.99) this, 
together with the formation of apparent horizons, provides strong evidence 
for the formation of a black string. 

As mentioned, for strong enough gauge fields, the previously described 
dynamics is severely affected. Figure 0] (left panel) shows the bubble area vs. 
proper time for different values of k. For large values he bubble "bounces" 
back and expands while for small ones the bubble collapses. There is a natural 
transition region separating these two possibilities. Tuning the value of k one 
can reveal an associated critical phenomena, the 'critical solution' being a 
member of the family of static solutions given by 
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Fig. 3. Rescaled Kretschmann invariant I 0N = Ir AH /12 vs. asymptotic time for 
m — 1.1 (solid line) and 1.99 (dashed line). The first non-zero values of the lines 
mark the formation of the apparent horizon. After some transient period, both lines 
approach the value of 1 suggesting a black string has formed. 



ds 2 = -V(r)dt 2 + IP-dr 2 + ^fldz 2 + r 2 V(r)dn 2 , (49) 



- F, v dx* A dx" = ±-V3r_(r + -r_)-^ , (50) 

where V(r) = 1 — r_/r and U(r) = 1 — r+/r. The parameters r_ and r+ 
(> r_) are related to the period of the z coordinate and to the ADM mass 
via P — 47rr + (l — r_/r + ) 3 / 2 and Madm = ^+/4. Since the quantities P 
and Maj^m are conserved, the member of the family of static solutions the 
dynamical solution approaches to can be determined a priori from the initial 
data. 

Figure^] (right panel) displays the time T defined as the length of asymp- 
totic time during which the bubble's area stays within 1% of the minimum 
value attained when the bubbles bounces back. This is a measure of how long 
the solution stays close to the static solution as a function of the parameter 
k. Empirically, we find the law 

T=-r+r\og\k-k c \+T 1 , (51) 

with a parameter r « 1.2 that does not seem to depend on the family 
of initial data chosen. This universality property is reinforced by the linear 
stability analysis of the critical solutions (|49I50|I performed in Ref. [35] where 
we prove that each solution has precisely one unstable linear mode growing 
like exp(/2t/r_|_) with a universal Lyapunov exponent of fl ~ 0.876. This 
explains the law (HU whh r = l/fi w 1.142. 
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Fig. 4. Left Panel. Area values vs. proper time at the bubble for different values 
of k and m = 1.1. By tuning the value of k appropriately, the amount of time 
that the area remains fairly constant can be extended for as long as desired. Right 
Panel The time T which is a measure of how long the solution stays close to the 
static solution vs. the logarithm of the difference between the parameter k and its 
critical value. A linear interpolation gives the value 7 « 1.2. 

3.2 Black holes 

As one of the applications that we have chosen to illustrate the use of the 
techniques previously discussed we consider here the evolution of single non- 
spinning black holes. Even when the data provided correspond to spherically 
symmetric and vacuum scenarios, as we will see, obtaining a long term stable 
implementation is not a trivial task. For additional information, and a more 
general treatment, we refer the reader to Ref. |33) . 

Formulation 

We adopt the symmetric hyperbolic family of formulations introduced in . 
This is a first order formulation whose evolved variables are given by 
{gij, Kij,dkij, a, A{\ with g^ the induced metric on surfaces at t = const, Kij 
the extrinsic curvature, dkij are first derivatives of the metric, dkij = dkgij, a 
is the lapse, and Ai are normalized first derivatives of the lapse, Ai — a~ 1 diCt. 

The Einstein equations written in this formulation are subject to the 
physical constraints, the Hamiltonian and momentum constraints, as well 
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as non-physical constraints, which arise from the variable definitions. The 
non-physical constraints are 



C Ai = Ai - N^diN = , C kij = d kij - d k9ij = , Q 



Ikij 



dud 



[ldk]ij 



= 0. 
(52) 

The constraints are added to the field equations and the spacetime constraint- 
functions {~f(t),((t),r](t),x(t),t;(t)} are introduced as multiplicative factors 
to the constraints. While these quantities are sometimes introduced as pa- 
rameters, we extend them to time-dependent functions. For simplicity in this 
work, we set £ = — 1. Requiring that the evolution system is symmetric hy- 
perbolic imposes algebraic conditions on these factors, and they are not all 
independent. If we require that all the characteristic speeds are "physical" 
(i.e. either normal to the spatial hypersurfaces or along the light cone), then 
we obtain two symmetric hyperbolic families. One family has a single free 
parameter, x(t)i 



Single constraint-function system < 



7 = 

c = 

v = 
£ = 



(53) 







and another symmetric system with two varying constraint-functions {n(t), 7(f) ^ 
-1/2}: 

rc= -1 



Two constraint-function system < 



7(2-??) 
l+2 7 



X 

2 

7^ -\ 



(54) 



Initial data and boundaries 

Initial data for a Schwarzschild black hole are given in In-going Eddington- 
Finkelstein coordinates. The shift ft 1 will be considered an a priori given 
vector field while the lapse is evolved to correspond to the time harmonic 
gauge with a given source function. This gauge source function is taken from 
the exact solution, such that in the high-resolution limit a = (1 + 2Af/V)~ 1 / 2 . 

Black hole excision is usually based on the assumption that an inner 
boundary (IB) can be placed on the domain such that information from this 
boundary does not enter the computational domain. This requirement places 
strenuous demands on cubical excision for a Schwarzschild black hole in Kerr- 
Schild, Painlevee-Gullstrand or the Martel-Poisson [HI] coordinates: the cube 
must be inside 0.37M in each direction. This forces one to excise very close 
to the singularity, where gradients in the solution can become very large, re- 
quiring very high resolution near the excision boundary to adequately resolve 
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the solution. This requirement follows directly from the physical properties 
of the Schwarzschild solution in these coordinates, and is independent of the 
particular formulation of the Einstein equations [Sj- 

With our current uniform Cartesian code, however, we do not have enough 
resolution to adequately resolve the Schwarzschild solution near the singu- 
larity. Thus, we place the inner boundary inside the event horizon, but out- 
side the region where all characteristics are out-going. The difference sten- 
cils are one-sided at the inner boundary, and no boundary conditions are 
explicitly applied. Testing various locations we find that placing the inner 
boundary at 1.1M gives reasonable results for the resolutions we are able to 
use, Ax = Ay = Az = M/5, M/10, M/20. We are working to resolve this 
inconsistency in our code by using coordinate systems that conform to the 
horizon's geometry. 

We performed numerical experiments with the outer boundary at three 
different locations, 5M, 10M and 15M. Boundary conditions for the outer 
boundary are applied using the orthogonal projection technique referenced 
above, by "freezing" the incoming characteristic modes. That is, their time 
derivative is set to zero through an orthogonal projection. This makes use 
of the fact that one knows that the continuum, exact solution is actually 
stationary. While this would not be useful in the general case, as we shall see, 
even in such a simplified case the constraint manifold seems to be unstable. 
We are currently working on extending the boundary treatment to allow for 
constraint-preserving boundary conditions and studying the well posedness 
of the associated initial-boundary value problem. 

Having set up consistent initial and boundary data, in a second order 
accurate implementation using the techniques mentioned in section II, we 
now concentrate on simulating a stationary black hole spacetime. As we will 
see below, even in this simple system, one encounters difficulties to evolve 
the system for long times. In particular, as has been illustrated in several 
occasions, the length of time during which a reliable numerical solution is 
obtained varies considerably depending on the values of the free parameters 
in the formulation. These parameters play no role at the constraint surface; 
however, off this constraint surface, these parameters have a sensible impact. 
Hence, at the numerical level -where generic data is only approximately at 
this surface-, it is necessary to adopt preferred values of these parameters. 
These, in turn, will depend not only on the physical situation under study 
but also in the details of the particular implementation (order of convergence, 
etc) . As we argued in Section II, the constraint minimization method provides 
a practical way to adopt these parameters. We next illustrate this in numerical 
simulations of Schwarzschild spacetime. 

Testing constraint minimization 

We concentrate here on black hole simulations performed using the symmetric 
hyperbolic formulation with two constraint functions. The single function 
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family and its disadvantages for constraint minimization are discussed in 

Ref. mm. 

Black hole numerical results 

As a first attempt to numerically integrate the Einstein equations, one could 
simply fix the parameters 77 and 7 to constant values. Lacking knowledge of 
preferred values for these parameters we might simply set 77 = and 7 = 0. 
Evolutions of the Schwarzschild spacetime for these parameter choices, how- 
ever, show that the solution is quickly corrupted, and the solution diverges. 
Figure El shows the error in the numerical solution with respect to the exact 
solution for three resolutions. While the code converges, the error at a single 
resolution grows without bound as a function of time. 



Convergence test, fixed constraint-functions ( y=0, r\=0) 
Errors with respect to the exact solution 




time [M] 

Fig. 5. Two-constraint-function family, with fixed values 7 = = 77, inner and 
outer boundaries at 1.1M and 5M, respectively. 



We now apply the constraint minimization technique to evolutions of a 
Schwarzschild black hole. The constraint functions rj(t) and j(t) will now vary 
in time, and both will be used to control the constraint growth. With two 
functions we can attempt to minimize changes in the functions themselves. 
This is advantageous because smoothly varying functions seem to yield better 
numerical results. Thus, rj(t) and j(t) are chosen at time step n+1 to minimize 
the quantity 

A := [7/(77 + 1) - 7)(ti)] 2 + [7(77 + 1) - 7 (n)] 2 (55) 

AT is nonlinear in 7 but linear in 77, allowing one to solve for 77 such that 
N = -aN, 
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(56) 

' X''(l + 27)+7Zx v > 

where, as in Section III, a is given by Eq. I|14|). 7 is chosen from some arbitrary, 
large interval. The corresponding rj given by Eq. i|5tj|) is computed, and the 
pair (77, 7) that minimizes A defined in Eq. (|55|l is chosen. 7 and 77 may 
be freely chosen, except that 7 7^ —1/2, giving two "branches": 7 always 
larger than -1/2, and 7 always smaller than -1/2. We have only explored the 
7 < —1/2 branch using the seed values T) = 0, 7 = — 1. In order to keep the 
variation of the parameters between two consecutive timesteps reasonably 
small, we have needed to set the tolerance value for the constraints energy 
roughly one order of magnitude larger than the initial discretization error, 
and n a to either 10 2 or 10 3 . This means that the constraints' energy, though 
in a longer timescale, will still grow. 



Dynamic minimization, boundaries at 5M 
Energy for the constraints and errors with respect to the exact solution 
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Fig. 6. This figure shows the constraint energy and the error with respect to the 
exact solution. Dynamic minimization is done with boundaries at 5M, Ax — M/5, 



10" 



and n a — 10 



The outer boundary is first placed at 5M. Figure shows the energy of 
the constraints and the error with respect to the exact solution. The corre- 
sponding constraint functions are shown in Figure [7| The large variation in 
the functions near the end of the run appears to be a consequence of other 
growing errors. In Figure |H1 the minimization is stopped at 750 M, and the 
functions are fixed to r; = —1.88, 7 = —1.00 for the remainder of the run. 
The solution diverges at approximately the same time. 

Another measure of the error in the solution is the mass of the apparent 
horizon, as shown in Figure |5| After some time, the mass approximately 
settles down to a value that is around 1.009M, which corresponds to an 
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Fig. 7. This figure shows the constraint functions for the run described in Figure|S| 
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Dynamic minimization, boundaries at 5M 
Keeping the constraint-functions fixed after t=750M 
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Fig. 8. Same as previous Figure, but keeping the constraint-functions constant 
after 750M. The figure compares the resulting energy for the constraints with that 
of the previous figure (shown at late times only, since because of the setup the runs 
are identical up to t = 750M). 



error of the order of one part in one thousand. For the higher resolution, 
the apparent horizon mass at late times becomes indistinguishable from 1M , 
given the expected level of discretization errors. 

The outer boundary is now placed at 15M. Figure ITU1 shows results for 
data equivalent to those discussed for Figure [S] The initial discretization 
value for the energy is 7.6459 x 10 -6 , and T = 10~ 5 , n a = 100 was used. The 
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Apparent horizon mass, boundaries at 5M 
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Fig. 9. This figure shows the black hole mass calculated from the apparent horizon 
with dynamic constraint-function values. The higher resolution simulation ran out 
of computing time. The apparent horizons were found using Thornburg's apparent 
horizon finder I3(i|. 



minimization of the constraint-functions is stopped at 450M, at which point 
the constraint- functions are approximately constant, and equal to 



Figure shows that the dependence of the lifetime on the location of 
the outer boundaries is not monotonic, as for this case the code runs for, 
roughly, lOOOAf , while with boundaries at 10M and 5M it ran for around 
700 M, and 800Af, respectively. A detailed analysis of such dependence would 
be computationally expensive and beyond the scope of this work, and may 
even depend on the details of the constraint minimization, such as the values 
for T and n a . However, comparing Figure 5 with Figures 6-9, we see that the 
constraint minimization considerably improves the lifetime of the simulation, 
as expected. 

4 Final Words 

We have chosen two problems to illustrate both the power of numerical sim- 
ulations of Einstein's equations and some of the difficulties encountered in 
obtaining accurate numerical solutions. This is especially relevant for black 
hole systems, where different poorly understood issues coupled to lack of suf- 
ficient computational power makes it much more difficult to advance at a 
sustained pace towards the final goal of producing a reliable description of 



j] = -1.35 x 10 



l 



7 



3.39. 



(57) 
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Boundaries at 1 5M 
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Fig. 10. Dynamic minimization done with boundaries at 15M, Ax = M/5, T = 
10 -5 , and n a = 10 2 . The constraint-functions are constant for t > 450M, where 
they are r\ = —0.135,7 = —3.389. Thus, the constraint functions do not respond 
when the code is about to crash. 



a binary black hole system. However, it is clear that goal outweighs these 
difficulties. As the bubble problem illustrates, a robust implementation was 
not only key to responding to open questions but also proved to be the way to 
observing other phenomena not previously considered. Not only did it show 
that a priori possible way to violate cosmic censorship is invalid, but it also 
revealed the existence of critical phenomena, which, in turn, can be used to 
shed further light in the stability of black string systems [32] ■ 

Fortunately, a substantial body of work in recent years has begun to 
address a number of these questions. A better understanding of the initial 
boundary-value problem in general relativity, advances in the definition of 



Recent analytical and numerical techniques 



27 



initial data and gauge choices coupled to several modern numerical techniques 
are having a direct impact in current numerical efforts. It seems reasonable to 
speculate that if this trend continues, the ultimate goal will be within reach 
in a not-too-distant future. 
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